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Abstract: 

Variational methods are used to calculate structural and thermodynamical properties of a titrating 
polyelectrolyte in a discrete representation. In the variational treatment, the Coulomb potentials 
are emulated by harmonic repulsive forces between all monomers; the force constants are used as 
variational parameters. The accuracy of the variational approach is tested against Monte Carlo 
data. Excellent agreement is obtained for the end-to-end separation and the apparent dissociation 
constant for the unscreened Coulomb chain. The short-range screened Coulomb potential is more 
difficult to handle variationally and its structural features are less well described, although the 
thermodynamic properties are predicted with the same accuracy as for the unscreened chain. The 
number of variational parameters is of the order of N 2 , where N is of the number of monomers, and 
the computational effort scales like N 3 . In addition, a simplified variational procedure with only 
two parameters is pursued, based on a rigid-rod approximation of the polymer. It gives surprisingly 
good accuracy for certain physical properties. 
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1 Introduction 



With increasing computer resources and refined algorithms it has now become possible to investigate 
structural and thermodynamic properties of polymer chains with several thousand monomeric units 
using simulation techniques. In particular, charged polymers have recently received an increased 
interest and a number of simulation studies have appeared in the literature [1, 2, 3, 4, 5, 6, 7, 8]. 
Most polyelectrolytes will under normal solution conditions exhibit an acid-base equilibrium, i.e. 
titratable groups in a polymer will exchange protons with the solution and the polymer net charge 
will vary as a function of the solution pH. This extra degree of freedom makes polyelectrolytes 
particularly versatile in many technical and biological applications [9, 10, 11, 12, 13, 14, 15]. In 
simulations, this phenomenon can be described by coupling the polyelectrolyte to an external proton 
bath of fixed chemical potential or pH; such grand canonical Monte Carlo simulations have been 
performed by several groups [16, 17, 18, 19]. 

Many polyelectrolytes undergo a dramatic structural change upon a successive ionization, i.e. the 
Coulombic repulsion between charged monomeric units forces the chain to adopt more extended 
conformations. This is particularly dramatic at low polymer and salt concentration. Sometimes the 
Coulombic repulsion is counteracted by attractive interactions leading to a helix-coil transition as 
a function of pH [20, 21, 22, 23]. In hydrophobically modified polyelectrolytes the attractive inter- 
actions tend to dominate and rather modest structural changes are seen upon ionization. Proteins 
are a special class of polyelectrolytes, which, although they may denature, still show rather small 
structural changes when going from the isoelectric point to either high or low pH. The competition 
between electrostatic monomer-monomer repulsion and specific attractive interactions has also been 
studied with simulation techniques [24, 17]. 

So far most polyelectrolyte studies have focussed on the behaviour of a single chain at infinite 
dilution. Additional salt has been introduced in an approximate way via screened Coulomb inter- 
actions. Finite polymer concentration and explicit salt particles impose considerable constraints 
on the chain lengths to be handled. A few studies, however, of non-titrating polyelectrolytes with 
explicit salt particles have been presented, but then only for fairly short chains [4, 25, 26, 27, 28]. 
Simulations with finite polymer concentration have also been carried through, but such studies are 
at present only feasible for a limited chain length [27, 28]. In this paper the accuracy of less computer 
demanding non-stocahstic approaches are investigated using Monte Carlo simulations. 

A number of approximate models for titrating polyelectrolytes have been suggested focussing on 
different phenomena and applying tools of varying mathematical sophistication. Existing theories, 
however, are seldom based on models that incorporate both a fully flexible chain and a discrete 
representation of the monomers. Here we suggest a variational approach, which has given excellent 
results for non-titrating polyelectrolytes [29, 30], based on a discrete representation of the polymer 
chain. That is, each monomeric unit is connected to its neighbours with a harmonic bond and may 
carry a unit charge depending on the solution pH. We apply the variational technique to a linear 
chain although the method is perfectly general and can be applied to any chain topology. The vari- 
ational solution to the titrating chain is obtained for different chain lengths and salt concentrations 
and the results are compared to results obtained via Monte Carlo simulations. 

The full variational solution is obtained at a much lower computational effort than what is required 
by simulations. As a further advantage it also provides the chain free energy. Still, the numerical 
procedure leading to the variational result is non-trivial, and we have also developed a simpler, but 
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of course less accurate, variational scheme based on only two variational parameters. 

This paper is organized as follows. In Sect. 2 we describe the details of the model we use for the 
titrating polymer. The variational technique of refs. [29, 30] is briefly reviewed in Sect. 3 together 
with some new numerical results for large N for non-titrating chains. In Sect. 4 the variational 
formalism is generalized to titrating chains, while the details of the Monte Carlo calculations can 
be found in Sect. 5. Sect. 6 contains the corresponding comparisons between Monte Carlo and 
variational results for chains in a salt-free environment as well as in salt solutions. Finally a brief 
summary and our conclusions can be found in Sect. 7. 



2 The Model 



The polyelectrolyte is regarded as an infinitely diluted polyacid in aqueous solution. The monomers 
form a linear chain where each monomer represents a titrating site that can be either protonated 
or deprotonated, i.e. be uncharged or carry one unit of negative charge. The chain is freely jointed, 
with neighbouring monomers connected by harmonic bonds. The implicit assumption is that the 
part of the underlying neutral polymer backbone that separates the neighbouring charged groups 
obeys a Gaussian distribution for its end-to-end separation. This is a reasonable approach for a 
chain where the separation of neighbouring charges is much larger than the persistence length of 
the underlying neutral chain. 

The solvent is treated as a dielectric continuum with a permittivity equal to that of water at room 
temperature. It also acts as a proton reservoir via a chemical potential given by 

ft= k B f \n 10 (pKo-pH) (1) 

where pH is that of the bulk, and pKo is the intrinsic pK a of a monomer. In the case of a salt solu- 
tion, an additional effect of the solvent is a Debye-Hiickel screening of the electrostatic interactions 
between all charged monomers, and the total interaction energy for N monomers becomes 

where x^- is the distance between monomer i and j, e is the electronic charge, e r is the dielectric 
constant of the solution (78.3 in all simulations), eo is the permittivity of vacuum and Zi the 
amount of charge on monomer i (either or —1). We use the tilde notation E, x^, etc. for physical 
quantities in conventional units, and reserve E, x^, etc. for dimensionless ones, which will be used 
in the variational formalism below. The force constant, k, is implicitly given through the input 
parameter fo = (e 2 /47re r eoA;) 1 / 3 , which is the equilibrium distance for a fully charged dimer and it 
is set to 6 A in all calculations. In the case of a 1 : —1 salt k is given by re = [2e 2 c s /e r eQNAkBT) 1 l 2 , 
where c s is the salt concentration (in mM), ks is the Boltzmann constant, Na is the Avogadro 
number and T is the temperature. 
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2.1 Dimensionless Formulation with Relative Coordinates 



Using dimensionless coordinates x^, given by = rox^, we may define a rescaled temperature T, 



T — 

kr 2 



and a likewise rescaled chemical potential /x, 



The negative exponent of the Boltzmann factor can then be written as 

EE 

~ = ™ ( 5 ) 

k B T T K ' 

with the rescaled energy given by 

E = 2 £ |x ili+ i| 2 + £ i x ..i + ME S * ( 6 ) 
where we have replaced Zi by more convenient {0, 1} variables Si = —Zi. 

In what follows, relative coordinates will mostly be used; instead of the absolute monomer positions 
Xi, the bond vectors r^, 

ri = x i+1 -Xi, i=l,...,N -1 (7) 

will be considered the fundamental variables. In this way complications due to the translational 
zero-modes are avoided; in addition, the convergence of the variational algorithm to be described 
below becomes considerably faster, especially at high temperatures. The energy of the chain then 
takes the form 



E(r, s) = E G + E C + E, = \ ^ r? + £ *' "^'^ + M E « 

1 = 1 a % 

where a runs over contiguous non-nil sub-chains, with 



corresponding to the distance vector between the endpoints, l a and r a , of the subchain. 



2.2 The Apparent Equilibrium Constant 

The average degree of dissociation, a = ]^X^( s i); can be interpreted in terms of an effective 
chemical equilibrium constant, pK 

P K=pH- lg (j^A (10) 
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Without the interaction, we would have 

=->(-£) 

A simple measure of the effect of the interaction on the chemical balance is then given by the change 
in pK: 

ApK = pK-pK = r- lg— — (12) 

F F F Tin 10 & l-a K ' 

3 Variational Treatment of Conformational Degrees of Free- 
dom 

3.1 Generic Formalism 

In this section we will focus on the conformational degrees of freedom, ignoring the titration; thus 
we consider a non-titrating polymer, with an energy obtained from eq. (8) by setting Si = 1 and 
disregarding the /x term: 

N — 1 

E(r) = E G + E c = -J2* 2 i+Y,— ( 13 ) 

? = 1 a 

A suitable variational approach to such a system is based on an effective energy Ansatz [31, 32, 29, 30] 

ij 

where define average bond vectors, around which Gaussian fluctuations are allowed, described 
by the symmetric, positive-definite correlation matrix Gij, the matrix inverse of which appears in 
the energy. 

Using this effective energy, the exact free energy F = —Tin Z of the polymer is approximated from 
above [33] by the variational one 

F = -TS V + (E) v > F (15) 

where Sy is the variational entropy, and (E)v is the average of the true energy in the trial Boltzmann 
distribution oc exp(— Ey /T). The parameters Gij and are to be determined so as to minimize the 
variational free energy F. The resulting effective Boltzmann distribution is then used to approximate 
expectation values (• • •) by effective (variational) ones (■■■)v- Thus, we have e.g. (rj)y = a^ and 
(ri ■ Vj)v = ■ &j + 3Gij. For potentials diverging like 1/r 3 or worse at short distances, (E)v will 
be divergent, and the approach breaks down. However, such potentials are not physical. 

At high T and/or small N, the resulting a^ will vanish. By setting a^ = in eq. (14), a restricted 
Ansatz is obtained, that in the screened case yields better results numerically than the a^ ^ case; 
this restricted version will be frequently used below. 

The minimization of F with respect to Gij and a^ gives rise to a set of matrix equations to be solved 
iteratively. These are considerably simplified, and the symmetry and positivity constraints on Gij 
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are automatic, if Gij is expressed as the product of a matrix and its transpose: 

JV-l 

/i=i 

The interpretation of the local parameter Zi is simple - it is a fluctuation amplitude for the ith bond 
vector r^. We can write 

r; = a; + ^ z ifl J ^ (17) 
/* 

where each component of £ 1Z 3 is an independent Gaussian noise variable of unit variance. 
Similarly, we have for a subchain 

r (7 = ^r i Ea (7 + ^ z all J M , (18) 

where a CT = X^go- a i anc ^ z f = Sign- z «- Thus, the noise amplitudes are additive. 
The matrix inverse of G can be similarly decomposed: 

G'ii = Wi • w,- (19) 
where Wi^ is the (transposed) matrix inverse of z^: 



Zi • Wj = 6ij (20) 



Note that Zi, Wi and z CT are vectors in R , not in R . 

In terms of and Zi, the variational free energy ignoring trivial additive constants (see ref. [30] for 
details) is given by 

P = -3Tlndetz + ^(3z l 2 +a l 2 ) (21) 

where 

= exp(a; 2 /2) erfc(a;/v / 2) (22) 
Setting a^ = for the unscreened case (re = 0) the variational free energy simplifies to 

F = -3T\ndetz+lj2 z l + ]fllL,j- ( 23 ) 

which very much resembles the energy of an (N — l)-dimensional Coulomb chain with bonds Zi, but 
with an extra entropy term (the first) preventing alignment of the ground state. 

The equations for a local extremum of F(a, z) are obtained by differentiation with respect to Zi and 

OF OF 

— = , — = 24 
ozi den 



5 



Due to the use of relative coordinates and of local noise amplitudes, a simple gradient descent 
method with a large step-size e can be used, that gives fast convergence to a solution of eqs. (24) 

8F 8F 
Azi = -e z —, Aa; = -e a — , (25) 
ozi den 

Further speed is gained by updating the reciprocal variables (arising from differentiating lndet z) 
using incremental matrix inversion (the Sherman-Morrison method [34, 30]) - the increment in Wj 
due to a change Azi is given by 

w,(w; • Az,) 

Aw, = - 1 (26) 

The resulting computational demand of the method is oc N 3 . 

The high and low T properties of the variational approach were analyzed in ref. [30]. In contrast 
to the case in MC simulations, the free energy is directly accessible with the variational method. 
Furthermore, the approach respects the virial identity, 2{Eq) — (Ec) = 3(iV — 1)T, in the sense that 
it holds also for the corresponding variational averages. 



3.2 Non-Titrating Variational Results Revisited 

The method, restricted to = 0, was extensively confronted with MC results in ref. [30]. Very 
good agreement was found for configurational quantities in the case of an unscreened Coulomb 
interaction (the error is well within the theoretical 11% limit [30]). In the screened case the method 
does not reproduce the MC results equally well although it gives a qualitatively correct picture of 
conformational properties. 

Prior to dealing with the titration case we will augment the comparisons in [30] on the end-to-end 
distance r ee with those arising from the a^ ^ solutions and also with those produced by the 
simplest possible variational Ansatz [35]. The latter is a highly constrained version of eq. (14), with 
Zi = and constant a^ = a, corresponding to a rigid rod. This simplification leads to simple scaling 
behaviour [35], 

r ee «tf(lntf) 1/3 (27) 
which seems to be approximately correct in MC calculations [30]; it is certainly correct at T — > 0. 
In fig. 1 we compare the results of the three different variational approaches (a^ = 0, a^ ^ 0, and 
the rigid rod) to MC data. The large N behaviour of the a^ = and a^ ^ variational curves 
is consistent with approaching the theoretical limits derived in ref. [29]; 11% and 0% respectively. 
The average monomer-monomer distances are also very well reproduced by both the a^ = and 
a^ ^ variational approaches. 

The simple rigid rod approach is more or less tailored to reproduce r ee , but it lacks the degrees of 
freedom to properly describe local properties along the chain. A similar approach was pursued in 
ref. [8] in which good results were reported, at temperatures significantly higher than those at which 
the calculations in this paper are performed. We expect the accuracy of the rigid rod approximation 
to increase with decreasing temperature. 

For screened Coulomb chains the a^ ^ results do not compare favourably with MC data - the 
chain tends to elongate more than the screened potentials call for. The same is true for the rigid 
rod approximation. 
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Figure 1: Relative difference (A=(Var-MC)/MC) in % for r ee between different variational versions 
and MC data for unscreened Coulomb chains as functions of N . The = 0, a^ ^ 0, and the rigid 
rod variants of the variational method are denoted by full, dashed and dotted lines respectively. 

4 Variational Approach to Titration 
4.1 Variational Treatment of Titratable Charges 

So far the variational approach has been confined to the situation with fixed identical charges along 
the chain. Next we generalize to allowing charge exchange with the solvent, with the total charge 
governed by a chemical potential. This amounts to considering Si in eq. (8) as dynamical variables, 
that each can be either or 1. Thus, suppressing for a moment the coordinate degrees of freedom, 
the system is isomorphic to an Ising spin system. 

E = -^^WijSiSj + [i^Si (28) 

ij i 

where the chemical potential /x has been introduced and coordinate dependencies etc. are lumped 
into the "couplings" Wij. Such systems have been subject to much attention in the solid state 
community, in particular for describing magnetic properties of so called spin glasses. A powerful 
alternative to computing thermodynamical properties of eq. (28) by means of MC techniques is the 
mean field (MF) approximation. This can be considered as a variational approach along the lines 
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above, with the variational energy Ansatz 

E v = -Tj2 u i s i (29) 

i 

where the coefficients Ui are the variational parameters. Minimizing the corresponding variational 
free energy one gets the MF equations 



Vi = g(ui) = ^(l + tanh (y)) 



, _ , , . . ,30 

2 V V 2 

dE 1 

Ui = (31) 

dvi T K ' 

which can be solved by iteration. The mean fields Vi have the interpretation of mean charges (si)v, 
while Ui are conventionally referred to as local fields. 

We will next merge the variational treatments of relative coordinates and charges treating the 
unscreened and screened cases separately. For simplicity we limit the presentation to the = 
solutions. The corresponding expressions for the ^ solutions can be found in Sect. 2. 



4.2 Unscreened Coulomb Chain 



For an unscreened Coulomb chain, the energy expression of eq. (8) simplifies to 

^E^ + E^ + mE* (32) 

where s\ a and s Ttr are the charges at the left and right endpoints of the sub-chain a, respectively. 
With a variational effective energy Ansatz 

Ev i T = \ E G ^ lri ■ r i - E uiSi ( 33 ) 

ij i 

including both coordinate and charge degrees of freedom, the variational free energy becomes, 
ignoring additive constants 

F = -3T In det z + T ^ { Ui Vi - In (1 + e Mi )} + | E z< ' + V ~ E + ^ E Vi ( 34 ) 

where Vi and Ui are related according to eq. (30); the first two terms give the entropy part —TSy, 
while the remaining terms correspond to (E)v • The equations for a local extremum of F(z, v) are 
obtained by differentiation with respect to Zi and Vi. One gets 

/2 x ^ vi v r z a , 
-E-4r^ = ° ( 35 ) 

<r3i " 

and 

M + Tui + J- V ^ = (36) 

respectively. 
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4.3 Screened Coulomb Chain 



The variational free energy for the screened case modifies to 

F = -3T\ndetz + T^2{ Ui Vi-\n(l + e Ui )} (37) 



i a \ a / i 



The derivatives then become 
OF 



<t3i 



1- k 2 zI) + k 3 zI-$(kz„)\ =0 (38) 



dF 12 1 \ 

— =M + Tui + ^«, ( Y~— - «*(««) = (39) 



4.4 Properties of the Variational Solution 
The Apparent Equilibrium Constant 

The variational approximation to a is simply given by 

a = JjJ2 Vi =¥= ff W ( 40 ) 
The local fields Ui can be expressed in terms of Vi (inverting eq. (30)) as 

Ui = g-\vi) =ln— ^- (41) 

1 - Vi 

Thus, the variational ApK is given by 

ApK\nlO = -^-g- 1 (^)) «-£-u (42) 
In the unscreened case (for simplicity), the last expression amounts to 

where the interpretation of ApK as an average energy cost per dissociated charge is very clear. 
The Virial Identity 

In the unscreened case, there is, also in the titrating version, a virial identity, 

2(E G ) - (E c ) = 3{N - 1)T (44) 
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obeyed by the exact thermodynamic ensemble. By taking the scalar product of eq. (35) with Zi 
and summing over i, we obtain 

- 3T(N - 1) + 3^z? - ^ = ( 45 ) 

and by a comparison with eq. (34), it is seen that the virial identity is indeed satisfied also for the 
variational expectation values. The same goes for the ^ solutions. 



Structure of the Effective Energy 



By taking the scalar product of eq. (35) with Wj, we obtain the effective force-constants 



™s'=««-sV;£ a ? ,L (46) 

while eq. (36) gives the effective local fields 

-7X = M+ J^£- (47) 

Thus, we can rewrite the variational energy, eq. (33), as 



(48) 



where the first and the last term reproduce the true bond and chemical potential terms of eq. (32), 
while the middle term emulates the effect of the Coulomb interaction on the charges (by appropriate 
local fields), and on the conformation (by suitable repulsive spring forces). 

The variational optimization thus forces the effective energy to have an intuitively appealing (and 
indeed very reasonable) structure. 



4.5 Algorithm Implementation 

For the numerical minimization of the free energy of eqs. (34) or (37) with respect to the variational 
parameters Zi (w^) and Vi (ui) a modified gradient descent is used; 

8F 8F 
Azi = -e z —, Ant = -e„ — . (49) 

OZi OVi 

where the use of the v (rather than u) derivative in the u-update implies a dynamical step size 

e,. 



v(l — v) 



(50) 



which is found to speed up the calculations. The Wi:s are obtained using incremental matrix 
inversion (eq. (26)). 

The complete algorithm looks as follows: 



10 



1. Initialize all Zi and Vi at random. 

2. Repeat until convergence: 
For all i: 

• Update Zi (if i < N) and Vi = g(ui) according to eq. (49). 

• Correct all Wj according to eq. (26) 

3. Extract Gij = Zi ■ Zj and compute variational averages of interest. 

Typical step-sizes are e z ps 0.15 and e u ps 0.5. The number of of computations for each iteration step 
is proportional to N 3 . An N = 80 system converges within about 100 iterations, which is somewhat 
slower than in the non-titrating case. 



4.6 A Simplified Variational Approach 

For a highly charged polyelectrolyte in the absence of salt, the conformation is strongly elongated. 
This fact motivates a simplified variational Ansatz, based on a rigid rod conformation and a constant 
degree of dissociation a along the chain. Such a simplified picture allows for analytical estimates of 
quantities of interest like ApK and r ee . To this end, consider the non-fluctuating limit Zi — > with 
a, = a = R/(-AT — 1) for all i, in which the variational Boltzmann distribution is given by 

exp{-E v /T) = exp (u£) S< ) II 6 ( r< " a ) ( 51 ) 

i 

Identifying a with g(u) gives for large N the variational free energy 

R 2 a 2 

F = NT {a ln(a) + (1 - a) ln(l - a)} + ttt; tt + ~^N{N - l)(ln N + 7 - 1) + N^a (52) 

z(iv — 1) it 

Here the sum 1/k has been approximated by In N+j, where 7 is the Euler constant. Minimizing 
F with respect to a and R gives for large N 

r ee = R « a 2/3 .ZV(m ^) 1/3 (53) 

(E c ) ^a 4/3 N{\nN) 2/3 (54) 

ApK « — — a^flnN) 2 ' 3 (55) 
F Tin 10 V ; K ' 

where /x has been eliminated in favour of a. This simplified model should give good results in the 
highly charged limit a — > 1, whereas its lack of fluctuations should give rise to poor performance in 
the a — > limit. 

A similar approximation can of course be used for a screened Coulomb potential. In this case, 
however, there is no analytical solution, but the relevant equations are readily solved by a simple 
numerical iteration scheme. 
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5 Monte Carlo Methods 



The Monte Carlo (MC) simulations were performed with the traditional Metropolis algorithm 
[36] in a semi-grand canonical ensemble. A single polyelectrolyte chain was simulated with the 
charges/protons moving between the monomers and an implicit bulk of fixed chemical potential. 

When a proton move is attempted, a monomer is picked at random and the charge state of the 
monomer is switched. The associated (free) energy change, AE, which determines if the move is to 
be accepted or rejected according to the Metropolis scheme, is the sum of the change AEc in the 
intramolecular Coulomb energy and a term ±/x that corresponds to the change in free energy for 
the acid-base reaction of an isolated monomer; the minus sign is used when the monomer is to be 
protonated, and the plus sign when it is to be deprotonated. 

In a single MC step a proton is only moved from the chain to the bulk or vice versa. Adding a step 
where a proton moves within the polymer does not affect the averages but increases the calculation 
time. 

When a conformational change is attempted, the associated energy change is given by the change in 
bond energy, AEq, plus the electrostatic energy change, AEc- The sampling of chain conformations 
is made highly efficient by using a pivot algorithm, which allows chain lengths of more than 2000 
monomers. The pivot algorithm was first described by Lai [37], and its efficiency for self-avoiding 
walks has been thoroughly discussed by Madras and Sokal [38]. 

A traditional move to update the coordinate variables is to attempt a translation of only one 
monomer at a time. The number of interactions that have to be calculated is of the order N 
for a highly charged chain and a large number of attempts per monomer is needed to generate 
independent chain conformations. In the pivot algorithm, however, each monomer i (except the 
first one) is translated in turn but together with the remaining semi-chain (monomers i + 1 to N). 
Furthermore, the semi-chain is then rotated as a rigid body around one of the coordinate axes with 
monomer i as origin. The number of interactions calculated in one step is of the order N 2 but 
independent conformations are obtained after only a few attempted moves, on the order of one per 
monomer or N in total. The net effect is a greatly reduced simulation time for a given degree of 
precision and a computational cost that grows approximately as N 3 . A completely different and 
even slightly more efficient procedure has recently been described by Irback [39]. 

A change in conformation is attempted once in every 20 steps; in the remaining steps, a change in 
the charge state is attempted. The total number of steps is around 10 8 . Every run is preceded by 
an equilibration of 10 5 -10 6 steps, where a change in conformation is attempted every other step, 
starting from a straight line with a charge on every other monomer. The simulations are faster 
at high jl, since fewer monomers are charged and thus the number of interactions that have to be 
calculated is smaller. 
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6 Results and Discussion 



6.1 Unscreened Coulomb Chain 

We now compare the variational and Monte Carlo results for r ee and ApK as a function of chain 
length and degree of dissociation. Fig. 2 shows an excellent agreement between the = variational 
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Figure 2: f ee /N as a function of a for for an unscreened polyelectrolyte chain. Monte Carlo results 
are denoted by crosses and the = and a^ ^ variational solutions are drawn as solid and dashed 
lines, respectively. The upper curves are for N = 1000 and the lower for N = 80 and the unit of 
length is A. 

results for r ee and MC data, with a maximal relative error of 8% for the N = 1000 chain at a = 1, 
well within the theoretical high- .AT limit of 11%. The difference between the variational and MC 
results increases with a as expected, but is indeed satisfactory for all cases studied. The a^ ^ 
variational solution is slightly inferior to the purely fluctuating solution, except for the N = 1000 
curve at high degree of ionization. Thus, we find a qualitative behaviour of the two variational 
solutions similar to that seen in fig. 1. 

The rigid rod approximation suggests that r ee /N should increase linearly with a 2 / 3 (ln N) 1 / 3 , which 
seems to be verified in fig. 3. Similar scaling behaviour for an analogous non-titrating polyelectrolyte 
has been derived in the literature [35, 7, 8]. We have, however, made numerical estimates of the 
exponents from logarithmic plots of r ee /N against a from both variational and MC data. For the 
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Figure 3: f ee /N as a function of a 2 / 3 (ln TV) 1 / 3 for different chain lengths (N = 20, 320, 1000) in the 
case of an unscreened Coulomb potential. Variational (a^ = 0) and MC results are denoted by solid 
lines and crosses, respectively. The unit of length is A. 



a dependence we find exponents in the range 0.80-0.85 from both approaches, which is significantly 
larger than the value of 2/3 predicted by the rigid rod approximation. The origin of the discrepancy 
is unclear. A similar discrepancy is seen for the In .AT exponent, which in the rigid rod case is given 
by 1/3, while numerical estimates from variational and MC data suggest a value of 0.6-0.7. One 
possible explanation could be that the chain expands via two different mechanisms. One is the 
expansion of each monomer bond, which should give rise to an In N exponent of 1/3 [30]; Another 
is the increasing alignment of monomer bonds with N , which reaches its limit when the coupling 
becomes strong and should level off for large N . 



Fig. 4 shows the shift in the apparent dissociation constant upon ionization of the chain and one 
finds that the = variational results differ significantly from the MC results for the longer chains. 
The a^ ^ results, on the other hand, are always in excellent agreement with the MC data with 
negligible differences for all systems studied. The largest error seen is of the order of one tenth of a 
pK unit. 



The global conformational properties of the polymer depends on the electrostatic coupling strength 
via an effective coupling, approximately given by 



a 4 / 3 N 



(56) 
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Figure 4: ApK as a function of a for different chain lengths with an unscreened Coulomb potential. 
Monte Carlo results are denoted by crosses and the = and ^ variational solutions are drawn 
as solid and dashed lines, respectively. The chain lengths are from top to bottom, N = 1000, 80, 20. 

The two variational solutions coincide for a weakly coupled chain, while at higher coupling two 
distinct solutions appear. This is seen in fig. 4 where the a^ ^ solution appears at a lower a value 
for the longer chains. Lowering the temperature (or decreasing ro) would have a similar effect. 

The asymptotic behaviour in the rigid rod approximation can be investigated by plotting ApK 
against a 1 / 3 (In N) 2 I 3 . Fig. 5 shows that the a dependence is not at all well described by the rigid 
rod over the parameter range studied and that it is only for very large N and a close to unity, 
that the slope seems to approach 1/3. A similar plot can be made for the In N dependence with 
a slightly better agreement and with a linear relation between hi ApK and ln(ln-AT). The slope 
is approximately 0.7 for a = 1.0, in good agreement with the rigid rod prediction of 2/3, but it 
increases with decreasing a. This is again consistent with two different expansion mechanisms. 

The rigid rod approximation gives an excellent description of the titration behaviour for a highly 
charged chain, while it slightly deviates from the MC results for short chains at low a. Fig. 6 shows 
that the rigid rod actually gives a better approximation to the MC data for ApK at high a values 
than the full a^ = variational solution. A similarly good agreement is found in fig. 7 where the 
MC and rigid rod numbers are virtually indistinguishable. The agreement between the rigid rod and 
the MC results holds for r ee and the apparent dissociation constant, but not for local properties like 
the monomer-monomer separation. The rigid rod does not distinguish between different positions 
along the chain and it also strongly underestimates the average monomer-monomer separation. For 
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In (a) 

Figure 5: In ApK as a function of In a for unscreened polyelectrolyte chains. From bottom to top: 
N = 20, 80, 320, and 1000. Variational (a^ = 0) and MC results are denoted by solid lines and 
crosses, respectively. 

such properties, the more sophisticated variational approaches are definitely superior and in general 
in good agreement with the exact MC results. 



6.2 Screened Coulomb Chain 

A discussion of the accuracy of the screened Coulomb approximation is beyond the scope of the 
present communication [4, 25]. Our principal aim is to investigate the variational techniques and we 
will use the screened Coulomb potential assuming that it contains the main physical features relating 
to screening. We will be primarily interested in r ee and ApK. These quantities will be functions of 
chain length and degree of ionization as before, but they will also depend on the screening parameter 
re. In a real solution re will contain contributions from added salt as well as any other charged 
molecule like the polyelectrolyte chain itself. Here we will neglect all these complications and we 
will refer to different screening conditions by stating the corresponding univalent salt concentration. 

A polyelectrolyte chain in a salt solution will be approximately independent of electrolyte concen- 
tration as long as the screening length is larger than the end-to-end separation and the screening 
will start to play a role when re -1 is smaller than or of the same order as r ee . 
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Figure 6: ApK as a function of a for an unscreened Coulomb chain. Monte Carlo results are denoted 
by crosses and the a^ = 0, a^ ^ and the rigid rod solutions are drawn as solid, dashed and dotted 
lines, respectively. The upper curves are for N = 1000 and the lower for N = 80. 

The accuracy of the variational result for r ee from the purely fluctuating solution (a^ = 0) dete- 
riorates with increasing chain length for a given salt concentration. The discrepancy to the MC 
data also becomes worse with increasing degree of ionization - see fig. 8. With increasing salt 
concentration, the accuracy will deteriorate up to some value, whereafter it improves. In the limit 
of very high screening the chain becomes perfectly Brownian, leading to an exact agreement. The 
variational predictions for r ee consistently overshoots; this can be attributed to an overestimate of 
the interaction with a Gaussian Boltzmann distribution. 

The variational solution with a^ ^ is absent at low a and high salt concentration. In the other limit 
there will always be two solutions and we find that the purely fluctuating one is the most accurate 
one in predicting the end-to-end separation. The a^ ^ solution on the other hand predicts much 
too large r ee - see fig. 8. The failure to describe the global structure of a screened chain is mainly due 
to the alignment of the bond vectors a^, which leads to the wrong asymptotic behaviour. The chain 
will expand like N instead of N 3 / 5 , which approximately seems to describe the MC results. The 
rigid rod approximation for the screened chain will of course show a behaviour for the end-to-end 
separation similar to that of the full a^ ^ solution. 

The apparent dissociation constant in the presence of salt (see fig. 9) is not well described by the 
fluctuating solution. Instead we find the a^ ^ solution, despite its obvious structural failure, to 
be superior and in fair agreement with MC data. Obviously the coupling between thermodynamic 
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Figure 7: ApK as function of N for a = 
results are denoted by crosses and the 
dashed and dotted lines, respectively. 



1.0 for an unscreened polyelectrolyte chain. Monte Carlo 
= 0, ^ and the rigid rod solutions are drawn as solid, 



derivatives, like the apparent dissociation constant, and the global structure is rather weak, and it 
is the local structure and fluctuations that determine the thermodynamics of a screened chain. The 
relative error is of the same size as for the unscreened chain - typically of the order of 10% or less. 
We find that ApK is consistently too large in the variational approximations, in accordance with 
the variational overestimate of the interaction. 

For the rigid rod approximation, the overall agreement with MC data is not as good as for the 
unscreened chain, and ApK is predicted to be essentially independent of a. However, as seen in fig. 
9, the numerical values obtained for a = 1 agree very well with simulated numbers, and the salt 
dependence of ApK for the fully ionized chain is accurately reproduced. 



7 Conclusions 



The variational calculations accurately reproduce structural and thermodynamic properties of a 
titrating polyelectrolyte interacting via an unscreened Coulomb potential. For a highly charged 
polyelectrolyte we find two distinct solutions to the Gaussian variational Ansatz, which at weak 
coupling collapse into one - a purely fluctuating solution with aj = 0. At strong effective coupling - 



18 



8 



6 



4 - 



O 




O 0.2 0.4 0.6 0.8 



Figure 8: f ee /N as a function of a for a screened Coulomb chain with N = 320. Monte Carlo results 
are denoted by crosses and the = and ^ solutions are drawn as solid and dashed lines, 
respectively. The salt concentrations are from top to bottom 0.001, 0.01 and 0.1 M and the unit of 
length is A. 

i.e. for large N , low T and large a - the a^ ^ solution always gives the lowest free energy, as well 
as the best approximation to structural and thermodynamic properties. At intermediate coupling 
strength, the purely fluctuating solution can produce superior structural data, while the apparent 
dissociation constant seems to be best described by the a^ ^ solution. 



With increasing electrostatic coupling, the chain becomes stiffer and more rodlike and a less extensive 
variational approach, like the rigid rod, becomes applicable. This can be solved analytically, and 
turns out to be fairly accurate both for thermodynamics and for certain structural properties. 



The screened Coulomb potential is more difficult to emulate with a Gaussian Ansatz, although 
reasonable predictions for the end-to-end separations are obtained from the purely fluctuating so- 
lutions. When the bond vectors a^ are non-zero they tend to align and the end-to-end separation 
increases linearly with N. For the purely fluctuating solution on the other hand, the end-to-end 
separation increases approximately as N ' 6 at high salt concentration, which is expected for a poly- 
mer with short range interactions. The apparent dissociation constant is well described also for a 
screened chain and we find that the a^ ^ solution, despite its structural shortcomings, is the best 
approximation for this property. 
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Figure 9: ApK as a function of a for a screened Coulomb chain with N = 320. Monte Carlo results 
are denoted by crosses and the = and ^ solutions are drawn as solid and dashed lines, 
respectively. The rigid rod results for a = 1 are indicated with arrows. The salt concentrations are 
from top to bottom 0.001, 0.01 and 0.1 M. 
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